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We study the combined effect of delayed Raman response and bit pattern randomness on pulse 
propagation in massive multichannel optical fiber communication systems. The propagation is de- 
scribed by a perturbed stochastic nonlinear Schrodinger equation, which takes into account changes 
in pulse amplitude and frequency as well as emission of continuous radiation. We perform extensive 
numerical simulations with the model, and analyze the dynamics of the frequency moments, the 
bit-error-rate, and the mutual distribution of amplitude and position. The results of our numerical 
simulations are in good agreement with theoretical predictions based on the adiabatic perturbation 
approach. 
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I. INTRODUCTION 



The interplay between noisy phenomena and nonlinear 
processes is a rich field of research that is of great inter- 
est in a variety of discipUnes including solid state physics 
[H, turbulence [1], and optics [1]. One of the most im- 
portant problems in this field concerns propagation of 
coherent patterns, such as solitons and solitary waves, 
in the presence of noise and/or disorder. An excellent 
example for systems where noise and nonlinear effects 
play an important role in the dynamics of coherent pat- 
terns is provided by fiber optics communication systems, 
which employ optical pulses to represent bits of informa- 
tion Q . It is by now well established that the parameters 
characterizing the pulses in fiber optics communication 
systems can exhibit non-Gaussian statistics [4, 5_, J, 0, [Sj . 
Yet, since optical fiber systems are only weakly nonlinear, 
it was commonly believed that the statistics of optical 
pulses is very different from the statistics encountered 
in strongly nonlinear systems, such as turbulence and 
chaotic flow, where intermittent dynamics exists. How- 
ever, a recent study of pulse propagation in optical fiber 
systems with multiple frequency channels in the pres- 
ence of delayed Raman response obtained results that 
stand in sharp contrast to this common belief [9|- This 
study focused on the interplay between Raman induced 
energy exchange in pulse collisions and randomness of 
pulse sequences in different frequency channels. Taking 
into account these two effects it was shown that the pulse 
parameters exhibit intermittent dynamic behavior in the 
sense that their normalized moments grow exponentially 
with propagation distance. Furthermore, it was shown 
that this intermittent dynamic behavior has important 
practical consequences, by leading to relatively large val- 
ues of the bit-error-rate (BER), which is the probability 
for an error at the output of the fiber line. 

The results of the study in Ref. Q were based on an 
adiabatic perturbation procedure that neglects radiation 



emission effects. However, these effects can be especially 
important for the fiber optics system under considera- 
tion. Indeed, as we shall see below, the interplay be- 
tween collision-induced energy exchange and randomness 
of pulse sequences can be described as an effective dis- 
order in the linear gain/loss coefficient, and the presence 
of gain can lead to instability with respect to emission of 
continuous radiation. Therefore, it is essential to obtain 
an improved description of pulse dynamics in the system 
that includes emission of continuous waves. In this pa- 
per we take this important task and derive a perturbed 
stochastic nonlinear Schrodinger (NLS) equation, which 
takes into account both changes in pulse parameters and 
radiation emission effects. We employ this model to an- 
alyze the dynamics of soliton parameters in comparison 
with the results of the simpler description of Ref. [9| and 
to draw conclusions on the possibility to observe inter- 
mittent dynamics in multichannel optical fiber commu- 
nication systems. The rest of the introduction is devoted 
to a summary of previous research on the effects of de- 
layed Raman response on soliton propagation in optical 
fibers. 

The main effect of Raman scattering on a single soli- 
ton propagating in the fiber is the self frequency shift. 
This effect, which is caused by energy transfer from 
higher frequency components of the pulse to its lower 
frequency components, was first observed experimentally 
by Mitschke and MoUenauer [l3| and explained theo- 
retically by Gordon 11]. Following this discovery, the 
impact of delayed Raman response on soliton propa- 
gation in optical fibers has drawn a lot of attention 
[l3, [H, [H, M, M, M, [H, Hi, HI • Most significantly, the 
influence on two-soliton collisions was studied by numer- 
ical simulations [ij, \1^ as well as by theoretical analysis 
[ij, ll3) 13 • These studies revealed that the main effect 
of a single two-soliton collision in the presence of delayed 
Raman response is an energy exchange between the col- 
liding pulses, which leads to a change of their amplitudes 



(Raman induced cross talk) [13|, [IJ, OJ, llll, lli ■ The fre- 
quencies of the two sohtons was also found to change as 
a result of the collision (Raman induced cross frequency 
shift) [ij, [l^, [ITI, liS] • Similar effects were recently stud- 
ied in collisions of ultra-short soliton pulses in photonic 
crystal fibers [2l|,|2i|. 

Raman induced energy exchange in pulse collisions can 
be beneficially employed in a variety of applications, in- 
cluding amplification in fiber lines [23, [2J| and in tunable 
laser sources [3|, [23] ■ However, it can also have nega- 
tive effects that impose severe limitations on the perfor- 
mance of multichannel communication systems. Indeed, 
it is known that the Raman-induced energy exchange in 
a single interchannel collision is independent of the fre- 
quency difference between the channels. Consequently, 
the magnitude of the induced energy shifts for a given 
pulse grows with the square of the number of channels 
[ISL l2a |. Thus, in a 100-channel system, for example, 
these effects can be larger by a factor of 2.5 x 10'^ com- 
pared with a two-channel system operating at the same 
bit rate per channel. Furthermore, since collisions with 
pulses from distant channels give the main contribution 
to energy shifts, a complete description of the dynam- 
ics must include interaction between pulses from all fre- 
quency channels. In contrast, effects of other nonlinear 
phenomena on pulse collisions are inversely proportional 
to some integer power of the frequency difference, and 
their cumulative influence can be adequately described 
by taking into account only a few neighboring channels 

Sua. 

Early studies of Raman cross talk in multichannel 
transmission systems focused on the dependence of the 
induced energy shifts on the total number of channels 
[281 [2g | . The combined effects of Raman cross talk and 
randomness of pulse sequences were also considered, and 
it was found that the probability distribution function 
(PDF) of pulse amplitudes is lognormal [H, H IM HI . 
However, these previous studies ignored several impor- 
tant properties of the system, which are essential for 
obtaining a correct dynamical model. First, all other 
nonlinear processes affecting pulse propagation, such as 
the Raman-induced self and cross frequency shifts were 
neglected. Second, strong coupling between amplitude 
dynamics and the dynamics of the other soliton param- 
eters, such as frequency, position and phase, was not 
taken into account. Consequently, only the amplitude 
PDF was calculated, whereas a correct evaluation of sys- 
tem performance requires calculation of the mutual PDF 
of the pulse amplitude and position. Third, most stud- 
ies considered only dynamic impact on performance of 
high frequency channels due to soliton decay, thus ig- 
noring potential negative consequences for intermediate 
and low frequency channels due to large position shifts 
induced by relatively large amplitude values. 

A more complete description of pulse propagation, 
which takes into account the three aforementioned fac- 
tors, was developed in Refs. [1, [33,|3J|- ^^ ^^^- O i* 
was shown that the coupling between frequency dynamics 



and amplitude dynamics leads to an exponential growth 
of the first two normalized moments of the self and cross 
frequency shifts with propagation distance. A perturbed 
NLS equation describing the combined effects of bit pat- 
tern randomness and Raman cross talk in a two-channel 
system was derived in Ref. [34]. Numerical simulations 
with the latter NLS model confirmed the analytic pre- 
dictions of Ref. [33] . Later on it was shown that the 
nth normalized moments of the self and cross frequency 
shifts increase exponentially with both propagation dis- 
tance and n^ ^. These results, combined with similar 
results for the normalized moments of the amplitude [33| | , 
imply that the soliton parameters exhibit intermittent 
dynamics in the sense that rare but violent events associ- 
ated with relatively large amplitudes and frequency shifts 
become important. Furthermore, it was shown that the 
dominant mechanism for error generation in the system 
at long propagation distances is related to the intermit- 
tent dynamic behavior and is due to the Raman-induced 
cross frequency shift ^ . In this process the error is gen- 
erated due to large values of the frequency and position 
shifts induced by large amplitude values. Thus, it is very 
different from the two mechanisms for error generation 
that are usually considered in fiber optics transmission, 
which are due to: (1) position shift with almost constant 
amplitude, (2) amplitude decay with almost constant po- 
sition shift. As mentioned above, the analysis in Ref. [3| 
ignored radiation emission effects, which can be impor- 
tant in massive multichannel transmission. In the current 
paper we take these effects into account and derive a per- 
turbed stochastic NLS model for propagation in a system 
with N > 2 frequency channels. We analyze the dynam- 
ics of the soliton parameters and the behavior of the BER 
by extensive numerical simulations with the model, and 
compare our results with the analytic calculations of Ref. 

The material in the rest of the paper is organized as 
follows. In Sec. [lT|A, we construct a perturbed stochas- 
tic NLS model describing soliton propagation in massive 
multichannel optical fiber transmission systems. The dy- 
namics of the soliton amplitude and frequency and of the 
BER are obtained in Sec. [IT| B by employing a standard 
adiabatic perturbation procedure. In Sec. IIIIl we analyze 
the results of numerical simulations with the perturbed 
NLS model and compare them with the predictions of 
the adiabatic perturbation theory. Section HVl is reserved 
for conclusions. 



II. A STOCHASTIC MODEL FOR PULSE 
PROPAGATION 

A. Derivation of the model 

Propagation of short pulses of light through an opti- 
cal fiber in the presence of delayed Raman response is 
described by the following perturbed NLS equation [3| : 



id,t^ + dfib + 2\'iP\^i^ = -e^VStlV'l', 



(1) 



where ip is proportional to the envelope of the electric 
field, z is propagation distance and t is time in the re- 
tarded reference frame. The term — e/j'09t|V'P accounts 
for the effects of delayed Raman response and en is the 
Raman coefficient 35). When tR — Q, the single-soliton 
solution of Eq. JT]) in a frequency channel (3 is given by 



%l;(i{t,z)=^-qp 



exp(tX;3) 

cosh(a:;^) ' 



(2) 



where xi^ = rjp {t-yp- 2(3z), xp ^ ap + [3(t - yp) + 

f 77^ — /3^ ) z, and rjp, ap and yp are the soliton amplitude, 

phase and position, respectively. 

Consider the effects of delayed Raman response on a 
single collision between two solitons from different fre- 
quency channels. For simplicity, one of the two channels 
is chosen as the reference channel with /3 = so that the 
frequency difference between the two channels is /3. We 
assume that e^ ^ 1 and 1/|/3| ^ 1, which is the typical 
situation in current multichannel transmission systems 
[36|. In addition, we assume that the two solitons are 
initially well-separated from each other in the temporal 
domain. Under these assumptions we can employ the 
perturbation proced ure, developed in Refs. [37, 38,^, 
and applied in Ref. [3J| for the case of delayed Raman 
response. Here we only give the outline of the calculation 
and refer the interested reader to Ref. [3J| for details. In 
accordance with this perturbative approach, we look for 
a two-pulse solution of Eq. ([T]) in the form 



Iptwo = ■00 + ■0/3 + 0, 



(3) 



where ^o and ipp are single-pulse solutions of Eq. ([T]) with 
< e/j <C 1 in channels and /3, respectively. The term 
on the right hand side of Eq. (I3|) is a small correction to 
the single-soliton solutions, which is solely due to collision 
effects. By analogy with the ideal collision case we take 
6 to be of the form 



= 00 + 0/3 + 



(4) 



where 0o and 0^ represent collision induced corrections 
in channels and /3, and the ellipsis represents higher 
order terms in other channels. Combining Eqs. ^ and 
(|4]) we see that the total pulse in the reference chan- 
nel is Vo°*"' = "00 + 00- We substitute the relations 
^ and ([U together with ■ipo{t,z) = ^o(a^o) exp(ixo), 
0o(i, z) = $o(a:^o)exp(ixo), ipp{t,z) = ■^p{xp)exp{ixp), 
and 0/3 (i, z) = ^p{xp) exp{ixp) into Eq. ([1]). The result- 
ing equation can be readily decomposed into an equation 
for the evolution of $0 and an equation for the evolution 
of ^p. We focus attention on $0 and remark that the 
calculation of $^3 is very similar. The equation for $0 is 
solved by integration with respect to z over the collision 
region. Carrying out this integration one obtains that 
the 0(ei?) effect of the collision on the reference channel 
soliton is given by 



A$ 



(1) 
01 ■ 



■r]psgn{P)eR'^o{xo), 



(5) 



where the first subscript in A$qj^ stands for the channel, 
the second subscript indicates the combined order with 
respect to both en and 1//?, and the superscript repre- 
sents the order in er. This OJ^r) effect corresponds to 
an amplitude change [H, [13, [13, llil 



A770 = 2rioripsgn{P)eH, 



(6) 



which is accompanied by emission of continuous radia- 
tion. In a similar manner, one finds that the effect of the 
collision in order e^//? is [34] 



A$ 



(i)_ 4z?7ffefl 



02 



dt'i' oixo). 



(7) 



A$Q2 corresponds to a collision induced frequency shift: 



A/3o = -{87^lvpeR)/{m), 



(8) 



which is also accompanied by emission of continuous ra- 
diation. 

Let us describe propagation of a reference channel soli- 
ton under many collisions with solitons from all other 
frequency channels in a system with 2N + 1 channels. 
We employ a mean- field approximation [34| , in which we 
assume that the amplitudes of the solitons in the other 
channels are constant. The random character of soliton 
sequences in different channels is taken into account by 
defining discrete random variables Qj , which describe the 
occupation state of the jth time slot in the «th channel: 
(^ij = 1 with probability s if the slot is occupied, and 
with probability 1 — s otherwise. It follows that the 
nth moment of Cij satisfies: (C") = s. We also assume 
that the occupation states of different time slots are un- 
correlated: (CyCi'j') = s^ if i 7^ i' and j ^ f . We de- 
note by A/3 the frequency difference between neighboring 
channels and by T the time slot width. Therefore, the 
distance traveled by the reference channel soliton while 
passing two successive time slots in the nearby channels 
is A.z^^' = T/{2Aj3). The 0{eR) effect of the collisions 
is taken into account by introducing a new perturbation 
term Si into Eq. ([1]). The term Si is obtained by sum- 
ming Eq. ([5]) over all collisions occurring in the interval 



,(1) 



Azc , and dividing the result by Azc 



(1) 



Si = ie_R*oe 



2X0 



^sg: 

i#0 



:n(/3.) J2 



j=(fc-i)i+i' 



.(!)' 



(9) 



where k — 1 and k are the indexes of the two successive 
time slots in the i = —1 channel, and the outside sum 
is from —N to TV. We decompose the disorder Qj into 



an average part and a fiuctuating part: ^i- 



c. 



where (Cy) = 0, {CijCi'j') = s{l- s)6u'Sjj', and Sw is the 
Kronecker delta function. Substituting Qj — s + Qj into 
Eq. ([9|) we obtain 



5^1 = 7^ > sgn(ft)H 



T 



JT^O 



-ieRaz)^o, (10) 



where the continuous disorder field ^(z) is 



a^) 



A, 



^ ki 



(11) 



i#0 



j = (fc-l)i+l 



Using Eq. pT|) and the properties of Cij one can show 
that (^(z)) = and (^(z)^(z')) = DnS{z ~ z'), where 
£)jv = N{N + 1)Z)2, -D2 = 2A/3s(l - s)r-\ and 5{z) is 
the Dirac delta function. Notice that the first term on 
the right hand side of Eq. (jlOp is zero due to symmetry. 
Even if this term is not zero, it can be compensated by 
appropriately adjusting the gain of the amplifiers. There- 
fore, the O(efl) effect of the collisions is described by 



^i = iej^^(z)^o- 



fl2) 



The 0{eji/(3) effect of the collisions is calculated in a 
similar manner. We first sum Eq. ([7]) over all collisions 
occurring within the interval Az, 



(1). 



1 



ki 



E 



c 



S2 = -cidt'i>o-4:€Rdt'i'oJ2 to ^^ 

^^0 l^'lj = (fe-l)»+l^^. 



IT, (13) 



where ci = (16A^e^s)/T. The second term on the 
right hand side of Eq. (fT3|) can be estimated as 

-8[Z?2i?w/(TA/3)]i/2ei?at*o, where Hn = Ef=i Vj- 
Consequently, for a typical multichannel system the co- 
efficient in front of e^dt'^o in this term is of order 1 or 
smaller, whereas for the first term this coefficient is of 
order TV. We therefore neglect the second term on the 
right hand side of Eq. ([T^ . and set S2 — —Cidt'^o- Using 
the fact that for a weakly perturbed soliton e^^°dt'i'Q = 
dt'ipo ~ iPoi^o we arrive at 



S2 = e'^"S2 = -cidti^o + ^Cl/3o^o, 



(14) 



where (3q is the frequency of the perturbed reference chan- 
nel soliton. Substituting Si and ^2 into Eq. ([!]) and 
replacing ip^ with ip we obtain 



id^i, + a,V + 2|vi V = -eflMIV'l' 



(15) 



which is the stochastic model describing propagation of 
the reference channel soliton in the fiber under many col- 
lisions. 



B. Statistics of soliton parameters and BER 
calculation 



The evolution of the parameters of the reference chan- 
nel soliton with propagation distance can be obtained 
by employing the standard adiabatic perturbation the- 
ory [40, [41] . Employing this perturbation procedure we 
obtain the following equations for the soliton amplitude 
and frequency: 



dm 
dz 



'^e-RS,{z)m{z) 



(16) 



and 



d/3o 8 4, . 



■^ciilliz). 



(17) 



Notice that the right hand side of Eq. (fT6|) is contributed 
solely by the second term on the right hand side of Eq. 
(flS)) . i.e., the term describing the Raman cross talk ef- 
fects. The first and second terms on the right hand side 
of Eq. (fTT]) describe the Raman induced self- and cross- 
frequency shifts, and are contributed by the first and 
third terms on the right hand side of Eq. (fT5|) . respec- 
tively. 

Integrating Eq. P^ over z we obtain 



m{z) = %(0) exp [2eBx[z)] 



(18) 



where x{z) = J^ dz'^(z') and r]o{0) is the initial ampli- 
tude. According to the central limit theorem, the PDF 
of x{z) approaches a Gaussian PDF with {x{z)) = and 
(a;^(z)) = Dnz. As a result, the PDF of the soliton am- 
plitude approaches a lognormal PDF: 

In^ [rjo/voiO]] 



F(7/o)-(8^Djve^z) 



-1/2 -1 

' Vo exp 



SDNe^fZ 



(19) 



The lognormal distribution is very different from the 
Gaussian distribution, and this difference is signifi- 
cant already in the main body of the distribution [34| . 
Moreover, the normalized moments of the lognormal 
PDF grow exponentially with propagation distance, from 
which it follows that the solijon amplitude exhibits inter- 
mittent dynamic behavior [3. 

The dynamic evolution of the soliton frequency is given 
by 



,(^) 



Poiz) = Pi;'iz)+p^'>iz) 



(c). 



where 



is the self frequency shift and 



P'o'iz) 



3T 



dz'4{z'), 



dz'r]l{z') 



(20) 



(21) 



(22) 



is) 



is the cross frequency shift. The nth moments of /Sq and 
/3o can be calculated from Q 



■0 (^)) 



n 



15 



e/?^o(0) 



dzm exp \?>2DNe\{2m - l)zm] , (23) 



and 



-\ciTll{Q) 



Wi'^^iz)) - 

X n / dz^exp [8DNe%{2m - l)z™] , (24) 
1 Jo 



where zq = z. By carrying out the integration in Eqs. 
(|23ll and ^ one can show that (/3o'^"(2)) and (/3^''^"(z)) 
are given by sums over exponential terms of the form 
if„exp [a^^^^^m^DArefjz], where a^") = 32, a^ = 8, 
< m < n, and the Km are constants. Furthermore, 
the leading contributions to the normalized moments 

(/?(^)"(^))/(/3(^n^))" and {Pi'^"iz))/{/3i''\z)r are expo- 
nentially increasing with both z and n?. As we shall see 
in the next Section, the normalized fourth moments of 
/3g , /3q and /3o increase much faster with increasing 
propagation distance compared with the normalized sec- 
ond and third moments, which is a consequence of the 
intermittent nature of the dynamics. 

In order to evaluate the system's BER we need to con- 
sider the main dynamical mechanisms leading to error 
generation. One mechanism, which has been widely stud- 
ied in relation with Raman cross talk ^, ^, ^, 3Q, _31, 
|32| . is due to pulse decay induced by loss of energy in 
collisions. This mechanism is associated with the small- 
T] tail of the amplitude PDF. Another mechanism, which 
has only recently been studied in relation with Raman 
cross talk, is due to the interplay between frequency (and 
position) dynamics and amplitude dynamics Q. In this 
case, the error is generated due to large values of the 
position shift, which are associated with large frequency 
shifts, and induced by relatively large values of the soli- 
ton amplitude. Notice that the lognormal statistics of 
the soliton amplitude leads to further enhancement of 
the BER contribution from the latter mechanism, since 
the large-?7 tail of the lognormal PDF lies above the cor- 
responding tail of the Gaussian PDF. We are therefore 
interested in the soliton position shift, which is given by 



yo{z) ^VoHz) 



-y^'\z), 



(25) 



where 



and 



y^oH^)--^-^ I dz'/ dz"r^tin (26) 



Jo 



dz' / dz"r]l{z') (27) 



yl^\z)--'-^ 



are the contributions from the self and cross fre- 
quency shifts, respectively. The position shift with 
a fixed amplitude rio{z) — 770(0) = 1 is yo{z) — 
yl'\z) + v\^\z), where y\,'\z) = -(8e,jz2)/15 and 
y'f'{z) = -{2,2N<iRSz'^)/(?,T). The relative position shift 
is Ayo(2) = ^yi"\z) -)- Ay^''(z), where ^y''^\z) = 
y'i\z) - y'^\z) and Ay(^)(^) = y'^^Xz) - y';^\z). We 
assume that yo can be compensated by employing filters. 
Therefore, the energy measured by the detector at a dis- 
tance z is 

/(770, A2/0) ^vl I dtcosh'^[T^o{t - A2/0)]. (28) 

J-T/2 



An occupied time slot is considered to be in error, if 
I{T]o,Ayo) < I{z = 0)/2 ~ 1. We estimate the BER by 
numerically integrating Eqs. ((26|l and ((27|) coupled to Eq. 
P^ for different realizations of the disorder (_{z) and cal- 
culating the fraction of errored occupied time slots. The 
^-dependence of the BER obtained by this calculation is 
described in Sec. IIIIl 



III. NUMERICAL SIMULATIONS 

In the previous Section we calculated the statistics of 
the soliton parameters and the BER by employing the 
adiabatic perturbation theory and neglecting effects as- 
sociated with emission of continuous radiation. We note 
that the latter effects can be particularly important for 
the system described by Eq. p3|) . Indeed, the second 
term on the right hand side of this equation has the form 
of disorder in the linear gain/loss coefficient. Such term 
can lead to instability with respect to emission of con- 
tinuous radiation, which is of second order in eji. It is 
therefore important to compare the results obtained in 
the previous Section by the reduced adiabatic method 
with results of numerical simulations with the more com- 
plete model, described by Eq. ([TS]). 

Notice that the fourth term on the right hand side of 
Eq. p5)) includes Po. Since both /3o and ci are of order er 
this term is of order e|j, whereas the other perturbation 
terms in the equation are of order e/j. Moreover, since (3o 
is a z-dependent random variable it is computationally 
complicated to solve Eq. (fTS]) in its exact form. To over- 
come this problem, we replace Po in Eq- (|15p with its 
value for the case where the amplitude is fixed and equal 
to 1: Poiz) = ~{32NeRSz)/{3T) - {8eRz)/15. Thus, the 
perturbed NLS which we solve numerically is 



+ieR£,{z)^ - 






(29) 



The initial condition is taken in the form of an ideal soli- 
ton: 2p{t,z = 0) = cosh~^(t), with ?7o(0) = 1, /3o(0) = 0, 
2/0(0) = 0, and ao(0) = 0. 

We perform Monte Carlo simulations with Eq. pS)) 
with about 5 x 10^ disorder realizations. The equation 
is integrated by employing the split-step method with 
periodic boundary conditions. Numerical errors resulting 
from radiation emission and the use of periodic boundary 
conditions are overcome by applying artificial damping 
at the vicinity of the boundaries of the computational 
domain. The size of the domain is taken to be —100 < 
t < 100 so that the absorbing layers do not affect the 
dynamics of the soliton pulses. The t-step and z-step are 
taken as At — 0.048 and Az = 0.001, respectively. 

We focus attention on a transmission system with 101 
channels operating at lOGbits/s per channel. It should 
be emphasized that state-of-the-art experiments with 
dispersion-managed solitons demonstrated multichannel 
transmission with 109 channels at 10 Gbits/s per chan- 
nel over a distance of 2 x 10'* km |43|. Several other 



experiments achieved total bit-rate capacities exceeding 
ITbits/s for shorter propagation distances [43, H, |45| . 
We use the foUowing set of parameters, which is similar 
to the one used in multichannel transmission experiments 
with conventional solitons ^36] , Assuming that T = b, 
A/3 = 10, s = 1/2 and 77i(0) = 1 for all channels, the pulse 
width is 20 ps, en — 3 x 10""', the channel spacing is 75 
GHz, and D2 — 1. Taking /32 — — Ips^/km, the soliton- 
peak-power is Pq = 1.25 mW. For these values the width 
of the lognormal PDF in Eq. (|19p . which represents the 



strength of disorder effects, is SD^ej^z = 1.8 x 10~ 2; 
for the reference channel. For z = 25, corresponding to 
transmission over 2 x 10* km, ^Dm^^z — 0.046. 

The z-dependences of the n = 2,3,4 normalized mo- 
ments of /3g , /3q , and /3o as obtained by numerical so- 
lution of the perturbed NLS are shown in Fig. [l]together 
with the results of the adiabatic perturbation theory. The 
numerical simulation results for (3^ and /3q were ob- 
tained by solving the reduced models 



0,V + 5,V^ + 21^-1 V = -eiiMIV-l' + Jefle(^)^, (30) 



and 



-cidtip + ici^o{z)'il! 



(31) 



with (3o{z) = -{32N€rsz)/{3T), respectively. The re- 
sults obtained by numerical solution of the perturbed 
NLS equation are in good agreement with those obtained 
by the adiabatic perturbation theory. Moreover, one can 

(s) (c) 

see that the fourth moments of /3o , /3o , and /3o increase 
much faster with increasing z compared with the second 
and third moments, in accordance with the intermittent 
nature of the dynamics. In addition, the normalized mo- 
ments of Pq grow faster than those of f3Q . This can 
be explained by noting that the rate of change of f]^ is 
proportional to tjq, whereas dfi^ /dz is proportional to 
?7o- Notice, however, that the values of the normalized 
moments of the total frequency shift j3q are very close to 
those of /3q . This is due to the fact that for the system 

described above /3q is typically much larger than /3q . 

The BER of the reference channel is calculated by the 
procedure described in Sec. II B. That is, we calculate 
the measured intensity using Eq. ([^5]) . and declare an 
occupied time slot to be in error, if I{'i]o,A.yo) < I{z — 
0)/2 ~ 1. The z-dependence of the BER obtained by 
numerical integration of Eq. (P9|) is shown in Fig. [2] 
together with the result obtained by employing the adi- 
abatic perturbation procedure. The agreement between 
the perturbed NLS simulations and the corresponding 
adiabatic theory calculations is good. Furthermore, it is 
seen that the BER attains relatively large values, which 
range from about 3 x 10"^ for z = 16 {X = 1.28 x 10" 
km) to about 10"^ at z = 25.0 {X = 2 x lO'' km). We 
remark that the BER values obtained by integrating Eq. 
(l5Tt , which takes into account only j/q , are very close to 




FIG. 1: Normalized moments of the reference channel soli- 
ton's cross frequency shift (a), self frequency shift (b), and 
total frequency shift (c) vs propagation distance z for a mul- 
tichannel system with 101 channels at 10 Gbits/s per channel. 
The solid, dashed-dotted, and dashed lines correspond to the 
n = 2, 3, 4 normalized moments obtained by the adiabatic 
perturbation method, using Eqs. (|22|l and (|24|l in (a), Eqs. 
(pT|) and ((23| in (b), and Eqs. (l20l)-((22l) in (c). The circles, 
squares and crosses represent the n = 2, 3, 4 normalized mo- 
ments obtained by numerical integration of Eq. (|3ip in (a), 
Eq. dSni in (b), and Eq. ^ in (c). 



the ones 

account 
between 
scale of 
with Eq 
and Eq. 
by notin, 



obtained by solving Eq. 

is) 



which takes into 

both tjq"' and yg''' . In fact, since the difference 
the two BER curves is indistinguishable on the 
Fig. [2l we choose to omit the result obtained 
, dSII)- The fact that the two models [Eq. ^ 
(I31|) ] give such close BER values is explained 
g that the cross frequency shift is typically much 
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FIG. 2: The z dependence of the BER for the reference 
channel in a 101-channel transmission system operating at 
lOGbits/s per channel. The circles represent the adiabatic 
perturbation prediction obtained by using Eqs. (|18|) . (|26|l . 
and H27I) . while the stars correspond to the result of numeri- 
cally integrating Eq. (I29|l . 



larger than the self frequency shift for the multichannel 
system considered here. 

As explained in Sect ion HIl and in Ref. Q, the two main 
mechanisms for error generation in the multichannel sys- 
tem are: (1) pulse decay, (2) large position shifts. While 
the first mechanism is associated with small pulse am- 
plitudes, the second one is predominantly due to large 
amplitudes. Hence, in order to better understand the 
roles of these two error-generating mechanisms one has 
to study the mutual PDF G{rio,Ayo). Figure [3] shows 
G{rio, Ai/o) for the reference channel soliton at the final 
propagation distance z = 25. The result obtained by nu- 
merical solution of Eq. (|29p (Fig. [3fa)) is in good agree- 
ment with the prediction of the adiabatic perturbation 
theory (Fig. [S^b)). Moreover, the mutual distribution 
function is strongly asymmetric about the Ayo = and 
770 = 1 axes. This asymmetry is a direct consequence of 
the strong coupling between position dynamics and am- 
plitude dynamics, as can be seen from Eqs. (pS]) and 
(P7)) . We emphasize that this behavior is very different 
from the one observed for soliton propagation in single- 
channel systems in the presence of amplifier noise, where 
the mutual PDF is approximately symmetric with re- 
spect to both Ayo — and % = 1 [6, _8|] . We also note 
that the mutual PDF shown in Fig. [3] is skewed toward 
larger 770 values, which can be explained by the skewed 
character of the lognormal distribution F{riQ). 



IV. CONCLUSIONS 

We investigated propagation of optical pulses in mas- 
sive multichannel optical fiber communication systems, 
taking into account the effects of delayed Raman response 
and bit pattern randomness. We derived a mean-field 



FIG. 3: The mutual PDF G(»7o, Ayo) for the reference channel 
soliton at 2: = 25 as obtained by numerical integration of 
Eq. (|29|l (a), and as predicted by the adiabatic perturbation 
method (b). 



description of the propagation, which is given by a per- 
turbed stochastic NLS equation that takes into account 
changes in pulse energy and momentum. This perturbed 
NLS model includes the effects of emission of continuous 
radiation, which were neglected in the simpler adiabatic 
perturbative approach used in Ref. '9|. Our numerical 
simulations show that the normalized moments of the 
soliton frequency shift grow exponentially with propaga- 
tion distance. Furthermore, the dynamics leads to rela- 
tively high values of the BER at intermediate and large 
propagation distances, and to an asymmetric form of the 
mutual PDF of amplitude and position. These results 
are in good agreement with results of the adiabatic per- 
turbation theory. We therefore conclude that the inter- 
play between Raman scattering and bit pattern random- 
ness plays a very important role in massive multichannel 
transmission systems. 

The exponential growth of the normalized moments 
of soliton parameters is indicative of intermittent dy- 
namics in the following sense: for certain realizations 
of pulse sequences (in the other frequency channels) the 
reference channel soliton can experience relatively large 
changes in its amplitude, which lead to relatively large 
position shifts and BER values. This dynamic behavior 
is quite surprising, since intermittency is usually associ- 
ated with strongly nonlinear systems such as turbulence 
and chaotic flow 2], whereas optical fiber systems are 
only weakly nonlinear. In a recent paper [4(?|, one of the 
authors suggested that this unexpected similarity might 
not be coincidental, but rather a consequence of the sim- 
ilarity between the dynamic behavior of the soliton am- 
plitude in the fiber optics system and the behavior of the 
local average of energy dissipation in turbulent fiow. 

We conclude by noting that the dynamic behavior de- 
scribed in the current paper is not limited to conventional 
optical solitons. Indeed, it is known that the effects of 
delayed Raman response on a single collision between two 
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dispersion-managed (DM) solitons are very similar to the 
effects in the conventional soliton case [1^. Therefore, 
we expect that similar results would hold for DM mul- 
tichannel transmission systems as well. A different type 
of nonlinearity that can lead to similar dynamics is due 
to nonlinear loss/gain. In this case pulse propagation 
is described by a perturbed NLS equation, in which the 
— CflVf^tl'^/'P term is replaced by ^pedV'I^V'j where tc is 
the cubic nonlinear loss/gain coefficient. It can be shown 
that the main effect of a fast collision in the presence 



of nonlinear loss/gain is a change in the soliton ampli- 
tude, which is given by an equation of the form ^ with 
sgn(/3)efl replaced by =F2ec/|/3|. If additional perturba- 
tions that affect the soliton frequency and position exist, 
the dynamics of the frequency or position will usually be 
coupled to the amplitude dynamics in a manner similar 
to the one described in Section II. Consequently, our re- 
sults should also be applicable for propagation of NLS 
solitons in systems with nonlinear loss or gain. 
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